#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin", "$FindBin::Bin");
use db_parameters;
use ortho_species;
use chromosomes;
use min_max;
use run_cmd;
use Data::Dump qw(dump);


use DBI;


print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_orthologous_exons


888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE


##########################################################################################

#	retrieve_gene_structure

#		retrieve data from ENSEMBL

##########################################################################################
sub retrieve_prot_id_to_cdna_hash($\%)
{
	my ($dbh, $prot_id_to_cdna_hash) = @_;
	print STDERR "\t\tRetrieve prot_id to cdna_hash...\n";
	my $data = $dbh->selectall_arrayref(<<"PL/SQLCMD") or die DBI::errstr;
	SELECT 
            prot_id,
            cdna_hash
        FROM 
            ens_id  NATURAL JOIN 
            cdna_sequences 
        WHERE 
            species = '$ortho_name1' or
            species = '$ortho_name2';
PL/SQLCMD
	for (@$data)
	{
		$prot_id_to_cdna_hash->{$_->[0]} = $_->[1];
	}
}
sub retrieve_ortholog_pairs($)
{
	my ($dbh) = @_;
	print STDERR "\t\tRetrieve ortholog pairs...\n";
	my $data = $dbh->selectall_arrayref(<<"PL/SQLCMD") or die DBI::errstr;
	SELECT 
            distinct ON (tc1.gene_id, tc2.gene_id) 
            tc1.gene_id, 
            tc2.gene_id 
        FROM 
            orthologs.ortholog_sets tc2 JOIN 
            orthologs.ortholog_sets tc1 USING (protocol_id, ortholog_id) 
        WHERE 
            tc1.protocol_id = $curr_protocol_id AND
            tc1.gene_id < tc2.gene_id AND 
            tc1.species <> tc2.species;
PL/SQLCMD
	return $data;
}



sub retrieve_coding_structure($$$$\%)
{
	my ($dbh, $species, $table, $columns, $prot_id_to_data) = @_;
	print STDERR "\t\tRetrieve coding structure from $table...\n";
	my $array = $dbh->selectall_arrayref(<<"PL/SQLCMD") or die DBI::errstr;
	SELECT
			prot_id,
			$columns
		FROM
			taxon.$table NATURAL JOIN 
            ens_id
        WHERE
            species = '$species'
		ORDER BY 
			prot_id, 
			rank;
PL/SQLCMD
	for my $line (@$array)
	{
		my ($prot_id, @data) = @$line;
		push(@{$prot_id_to_data->{$prot_id}},  [@data]);
	}
}

sub retrieve_chromosome_strand_by_prot_id($$\%)
{
	my ($dbh, $species, $data) = @_;
	print STDERR "\t\tRetrieve chromosome and strand...\n";
	my $array = $dbh->selectall_arrayref(<<"PL/SQLCMD") or die DBI::errstr;
    SELECT 
            DISTINCT ON (prot_id) 
   			prot_id,
            'chr' || chromosome,
            strand
        FROM 
            ens_id 
            NATURAL JOIN taxon.gene_loci
        WHERE
            species = '$species';
PL/SQLCMD
	for my $line (@$array)
	{
		$data->{$line->[0]} = $line->[1]."\t".($line->[2] ? '+': '-');
	}
}


sub retrieve_gene_structure($$$\%)
{
	my ($dbh, $species, $filename, $prot_id_to_cdna_hash) = @_;

	open(FILE, ">$filename") or die;

	
	my (%prot_id_exons, %prot_id_introns, %prot_id_chromosome_strand);
	retrieve_coding_structure($dbh, $species, 'prot_coding_exons', 'length, frame, exon_id', %prot_id_exons);
	retrieve_coding_structure($dbh, $species, 'prot_coding_introns', 'length', %prot_id_introns);
	retrieve_chromosome_strand_by_prot_id($dbh, $species, %prot_id_chromosome_strand);

	my @ret_array;
	for my $prot_id (keys %prot_id_exons)
	{
		my ($chrm_strand) = $prot_id_chromosome_strand{$prot_id};
		die $prot_id unless defined $chrm_strand;
		my $exons = $prot_id_exons{$prot_id};
		my $introns = $prot_id_introns{$prot_id};
		$introns = [] unless (defined $introns);

		# handle first exon
		die $prot_id unless @$exons;
		die $prot_id unless defined $exons->[0][0];

		# merge "NONE" into first exon
		if ($exons->[0][2] eq 'NONE')
		{
			$exons->[1][0] += $exons->[0][0];
			shift(@$exons);
		}
		# merge "NONE" into last exon
		if ($exons->[-1][2] eq 'NONE')
		{
			$exons->[-2][0] += $exons->[-1][0];
			pop(@$exons);
		}

		
		my $cumm_exon_pos = $exons->[0][0];
		my $cumm_pos = $cumm_exon_pos;
		#	prot_id chrm strand frame rank cumm.start cumm.start + len cumm.intron_start + len
		print FILE "$prot_id_to_cdna_hash->{$prot_id}\t$chrm_strand\t0\t1\t0\t$cumm_exon_pos\t0\t$cumm_exon_pos\n";

		# get rid of first exon now we have printed it out
		shift(@$exons);

		unless (@$exons == @$introns)
		{
			print STDERR (dump($introns))."\n";
			print STDERR (dump ($exons))."\n" ;
			print STDERR scalar (@$exons). ",". scalar (@$introns). "\n";
			die $prot_id;
		}

		# remainding exons
		for my $i (0.. (@$exons - 1))
		{
			my ($exon_len,  $frame) = @{$exons->[$i]};
			my ($intron_len) = @{$introns->[$i]};
			# this will blow up if prot_coding_introns table has null length
			# in a previous bug in 23_table_prot_coding_exon/table_prot_coding_exons.pl
			#unless (defined $intron_len)
			#{
			#	print STDERR (dump($introns))."\n";
			#	print STDERR (dump ($exons))."\n" ;
			#	print STDERR scalar (@$exons). ",". scalar (@$introns). "\n";
			#	die $prot_id;
			#}
			$cumm_pos += $intron_len;
			print FILE join ("\t",  $prot_id_to_cdna_hash->{$prot_id}, 
									$chrm_strand,
									$frame,
									($i+2),
									$cumm_exon_pos,
									$cumm_exon_pos + $exon_len,
									$cumm_pos,
									$cumm_pos + $exon_len), "\n";
			$cumm_exon_pos += $exon_len;
			$cumm_pos += $exon_len;
		}
	}
}

sub retrieve_gene_to_transcript_map($$)
{
	my ($dbh, $species) = @_;


	#
	#	ENS table name
	#
	my $array = $dbh->selectall_arrayref(<<"PL/SQLCMD") or die DBI::errstr;
    SELECT
            gene_id, 
            start - 1, 
            finish, 
            strand, 
            'chr' || chromosome,
            prot_id
        FROM
            ens_id NATURAL JOIN 
            taxon.gene_loci
        WHERE
            species = '$species'
PL/SQLCMD
        return $array;
}
#sub retrieve_gene_to_transcript_map($$)
#{
#	my ($mysql_dbh, $species_ens_abbrev) = @_;
#
#
#	#
#	#	ENS table name
#	#
#	my $table_name = $species_ens_abbrev . '_gene_ensembl__transcript__main';
#	my $array = $mysql_dbh->selectall_arrayref(<<"PL/SQLCMD");
#	SELECT
#			gene_stable_id,
#			gene_chrom_start -1,  
#			gene_chrom_end,
#			chrom_strand,
#			CONCAT('chr', chr_name), 
#			translation_stable_id
#		FROM
#			$table_name
#		ORDER BY 
#			gene_chrom_start,
#            gene_stable_id,
#			translation_stable_id
#PL/SQLCMD
#        return $array;
#}



	#
	#	Get gene to transcript map data
	#
	my $dbh = connect_to_panda();
	$dbh->{RaiseError} = 0;
	use get_psql_time;
	my $start_time = get_psql_time($dbh);

	my %prot_id_to_cdna_hash;
	retrieve_prot_id_to_cdna_hash($dbh, %prot_id_to_cdna_hash);

	
	my $species_index = 0;
	my @species_data = ([$ortho_name1, $ortho_prot1], [$ortho_name2, $ortho_prot2]);
	for my $species_datum(@species_data)
	{
		++$species_index;
		my ($species,  $regex) = @$species_datum;

		print STDERR "\t$species:\n";
		#
		#	prepare gene structure
		#
		{
			print STDERR "\tRetrieving gene structure for $species...\n";
			retrieve_gene_structure($dbh, $species, 
							"$FindBin::Bin/gene_structure_".$species_index.".cds", 
							%prot_id_to_cdna_hash);
		}

		#
		#	prepare gene to transcript map
		#
		{
			print STDERR "\tRetrieving gene to transcript map for $species...\n";
			my $data = retrieve_gene_to_transcript_map($dbh, $species);
			open(FILE, ">$FindBin::Bin/genes2transcripts_".$species_index.".map") or die;
			#print FILE "stable_id\tseq_region_start\tseq_region_end\tseq_region_strand\tchrname\tstable_id\n";
			for my $line (@$data)
			{
				my $prot_id = $line->[-1];
				my $cdna_hash = $prot_id_to_cdna_hash{$prot_id};
				die $prot_id unless exists $prot_id_to_cdna_hash{$prot_id};
				print FILE join ("\t",  @{$line}[0..4], $cdna_hash), "\n";
			}
		}

		
	}
	{
		open ORTHO_PAIRS, ">$FindBin::Bin/all.orthologs.gene_pairs" or die;
		my $ortholog_pairs =  retrieve_ortholog_pairs($dbh);
		for (@$ortholog_pairs)
		{
			print ORTHO_PAIRS $_->[0]."\t".$_->[1]."\n";
		}
		close ORTHO_PAIRS;
	}


	# change to makefile directory
	use Cwd qw(abs_path getcwd); # 'abs_path';
	my $curr_dir = getcwd();
	(my $my_dir = $FindBin::Bin) =~ s/(.*)\/$/$1/;
	#print STDERR $my_dir, "\n";
	chdir($my_dir);
	
	system ("gmake -f $FindBin::Bin/makefile -j 100 ".
	 		"PROT_SEQ1=$FindBin::Bin/../sequences/$ortho_name1/prot.mapped.fa ".
			"CDNA_SEQ1=$FindBin::Bin/../sequences/$ortho_name1/cdna.mapped.fa ".
			"PROT_SEQ2=$FindBin::Bin/../sequences/$ortho_name2/prot.mapped.fa ".
			"CDNA_SEQ2=$FindBin::Bin/../sequences/$ortho_name2/cdna.mapped.fa "
			);

unlink("$FindBin::Bin/gene_structure_1.cds");
unlink("$FindBin::Bin/gene_structure_2.cds");
unlink("$FindBin::Bin/genes2transcripts_1.map");
unlink("$FindBin::Bin/genes2transcripts_2.map");
unlink("$FindBin::Bin/all.orthologs.gene_pairs");

print STDERR "\tCompleted\n";
# restore directory
chdir ($curr_dir);

log_pipeline_stage($dbh, $curr_protocol_id, 38, 'table orthologous_exons', $start_time);
$dbh->disconnect();
